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THEORETICAL STUDIES OF A MOLECULAR BEAM GENERATOR 

MOLECULAR BEAM GENERATOR MODEL 

The following is a proposed baseline model that is being develope for the simulation of 
hydrodynamic generator, which can be converted at a later date to a magnetohydrodynamic 
MHD thruster by adding the necessary electric and magnetic fields. The following development 
will include the electric and magnetic terms, however, the initial computer program will not 
include these terms. The analysis that follows is for a one species, single temperature model 
constructed over the domain D defined by the region enclosed by ABCDEF illustrated in the 

figure 1. 



Figure 1. Geometry of thruster MHD model. 


CONTINUITY EQUATION 


The continuity equation expresses conservation of mass and is given by 




(i) 


i 


where p = p(r, 0, z,t) is the density of the gas, and V = V r e r + Vg eg + V z e z is the velocity. In 
cylindrical coordinates the equation (1) has the form 

dp I cHspVr) I djpVg) djpVz) _ q ( 2 ) 

~5t + r dr r d9 dz 


CONSERVATION OF MOMENTUM 


The equation for conservation of linear momentum is given by 

+ P(V^)V = j^F, 
m ,= 1 


( 3 ) 


where represents a summation of body forces per unit volume acting upon a control 

volume within the domain D. We consider initially the pressure force 


F\ = -VP ( 4 ) 

where pressure and density are related by the equation of state gas law P = p* RT , where p* is 
the density in mole/m 3 , i.e. p*W = p where W is the molecular weight in kg/mole. The force 

due to viscosity is 


F 2 = ti { V 2 u + V(V • V)} - ^ V(T ? V • V) + 2 (Vrj ■ V)V + Vrj x (V x V) 


( 5 ) 


where 


,0 5 (Mmk B T\ l/2 f2k B T\ 2 nrrb ,2 

^ = 1 - 2510 ' s^a{— ) {-?-) =ClT 


(6) 


is the plasma viscosity, with 


A = 2(log(l + ft 2 ) 


ft 


r) 


1 + o; 2 

a constant which depends upon the ionization factor ft, (ft = 1 for a fully ionized gas). The 
additional constants are M the ion mass (Kg), m the electron mass, e the electron charge, k B 
Boltzmann’s constant, T is the absolute temperature, and n; = 1 for a singly ionized plasma. 
For ft = 0, we employ an empirical curve fit for the viscosity as a function of temperature. The 

magnetic force is given by 

p 3 = / x B ( 7 ) 


2 


2 


where J is the current density. The gravitational force is given by 


The electric force is given by 

II 

< 

All additional forces are represented by 

F b = eE. 

n 

and are neglected for the present. 

CONSERVATION OF ENERGY 

E* 

i =6 


Representing the internal energy by u = C p T where C p is the specific heat at constant pressure 
and T is the absolute temperature, the energy equation can be written in the form of an energy 
balance as n 

d(C£) + ^ . V )(C P T) = V(K t VT) + *i (8) 

& i-\ 

where both the specific heat C p and thermal conductivity Kt are treated as functions of 
temperature T. The thermal conductivity K T of the medium is given by the Spitzer-Harris 

relation 10^/9 

_ 4 ,4io-.. T V2 ^ (9) 


™ , r 1.22 10 3 n^' 

23 - log ^^72 


and n is the plasma number density in particles/m 3 . In addition to the heat loss term the right 
hand side of equation (8) contains the terms 

fa =(E + V x B)- f-rjV ■ E ( 10 ) 


which represents joule heating, 
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which represents viscous dissipation. Here rj and A are viscosity coefficients satisfying \ + jt) = 0. 
In addition there is the radiation loss term. Various forms of the radiation term exists in the 
literature. As a first approximatrion we take the radiation loss term from reference 11 which can 

be expressed ^ 

fa = VAfiVt-aT 4 ) 

where \ R is the Rosseland mean free path (\ R = l/a R where a R is the Rosseland absorption 
coefficient (cm -1 )), and cr is the Stephan Boltzman constant. The remaining terms EtU 
represents additional energy considerations which are initially neglected. 

MAXWELL’S EQUATIONS 

Maxwell’s equations in the MRS Rational system of units can be expressed 


Gauss’s law for magnetism 

VB = 0 

(12) 

Gauss’s law for electricty (Coulomb’s law) 

V ■ D = pe 

(13) 

Ampere’s law 

dD 

Vx ' ff = J + -a 

(14) 

Faraday’s law 

- SB 

V * E= Sf 

(15) 


CONSTITUTIVE EQUATIONS 

Assuming an isotropic, homogeneous medium we adopt the constitutive equations 

D = e E and B = fi H. 


OHM’S LAW 

Ohm’s law is written in the form 

J = a(E + VxB)-^(JxB) (16) 

where J is the current density, E is the electric field, B is the induced magnetic field and a is 
the electric conductivity with units of mho/m and fi is the Hall parameter given by (reference 1) 

fi = 9.6(10 1 6)(T 3/2 R/Zn log A) 


4 


4 


with Coulomb logarithm given by 


log A w 23 - log(1.22 x 10 3 n 1 ! 2 /T*! 2 ) 


with n the plasma number density. 


ELECTROMAGNETIC FIELD EQUATIONS 


Neglecting the displacement current modifies Ampere’s law to 

VxS = fxJ. 


(17) 


Assumming that the charge density is constant implies the equation of continuity of charge 
is V • J = 0. (Note that the divergence of equation (17) also gives this result.) From Jackson, 
reference 10, along with neglecting the displacement current, it is appropriate to ignor Coulomb s 
law as its effects are negligible. We thus obtain the electromagnetic field equations 


VxB = /iJ 


Vx£ = 


dB_ 

di 


(18) 


Using equation (17) in Ohm’s law we solve for E and write 

B = aVxB-yxB + /?(VxB)xB 


(19) 


where /3 = 0,/iJ.Ba and a = l/api. We substitute the results from equation (19) into Faraday s 
law and write 

= V x (aV xB)-Vx(V xB) + Vx £(V x B) x B (20) 

dt 

and since /? is a function of T we find 

- v x (aV x B) - V x (U x B) + V/? x (V x B) x B + /3V x (V x B) x B (21) 
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SUMMARY OF BASIC EQUATIONS USED FOR MODELING 


Continuity 

! + v (,vo = ° 

Momentum 

^ +/>(?■ v) ? = E £ 

m i = 1 

Energy 

n 

p^P- + p(V • V)(C P T) = V(K t VT) + Y, 
at ,=i 


Electromagnetic field equations 


dB 

dt 


= V x (ftV x B) - V x (V x B) + V x 0(V x B) x B 


This produces a system of eight simultaneous partial differential equations m the eight 
unknowns 

Br,B e ,B z ,p,V r ,V 0 ,V z ,T. 

Throughout the calculations the following quantities can be generated in terms of the above 
variables. ^ 

U ~p VX§ _ ( 22 ) 

E = apJ - V x B + J x B) 


SCALAR FORM OF FIELD EQUATIONS 

Assu min g symmetry with respect to the 6 variable, all derivatives with respect to 6 are 
neglected. The following set of scalar equations then results 

Continuity 

dp l djrpVr) d[pV z ) 

dt r dr dz 


6 
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Momentum 
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Electromagnetic field equations 

dBr d 2 Br . d 2 B z . . 


= -“lP“ + “a ^7 + ’ 4 "ar + r a* y " a* 


„ f„ a 2 fl„ , dB » dB e , u f a 2 B« 1 aB«\ 

~ /) [ B *"a?' + - a 7 "ar + r Va--az + r a* ; + a* 


dz V dr 


a/? r D as* 


_ * dz 

a 2 B fl a 2 B fl a dBe a n v as 2 „ ^ 


aBg Be 
dr r 


dB tt „ av 2 Tr aB* . _ av r aB r av* 

+ + B# a7 + Vr lT + B * aT - v, -aT " ^ 

[ /d 2 B r d 2 B 2 \ dB z (dB r dB z \ 2B 0 dBe 

+ /? M"dP~ drdz + dz \ dz dr ) r dz 


z \ dz 2 drdz 
d 2 B r _ d 2 B 2 \ 
drdz dr 2 ) 

' /dB r ae, 
* V dz dr 
d 2 By d 2 B 2 


dr V dz 


B# "az +Br 


dB r l dBr as/\l , d0 f„ dB„ , „ ( dB t dB, 

dz dr 


-W - B ° 


dBe Be_ 
dr r 


d 2 B z | a dB r a dB 
dr 2 r 3z r dr 


' dz L V az ar ) " \dr 

dB z _ d 2 B r d 2 B z c*dB r Q 

~df ~ a drdz a dr 2 r dz r 

dB z „ dV r 1 /T , D T/ D s 

+ v r ^ + b 2 + - (v r B 2 - v 2 B r ) 

dr dr r 

. „r„ a 2 B« . aft as, , „ fa 2 ^ 


£ -^- B 


, au 

r 5r 


aw i 


, a n Z-Z± A 0 + B — — Z. + -— 

+ ^ ^ drdz + dr dz r \ dr 2 r dr 


d 2 B e 2dB 0 \ , dB r /dB* B e \ B 2 dBe 


dr V dr 


B z ^Z + Br 

dz 


dBe B# 

dr r 


These equations are subject to certain boundary and initial conditions which are now 
discussed. 

BOUNDARY AND INITIAL CONDITIONS 

With reference to the figure 1, the line AF has the input conditions 

p = Pq — constant 
T = To = constant 

V r = V e = 0 

V z = Vq — constant 
dBr _ dBe_ _ dBz 
dz dz dz 
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Due to symmetry considerations the line AB has the center line boundary conditions 


dp 

dr 

dT 

dr 

dVr 

dr 

dVe_ 

dr 

dVz_ 

dr 


0 

0 

0 

0 

0 


dBr 

dr 

dBe 

dr 

dB z 

dr 


= 0 
= 0 
= 0 


The far field conditions along the line BC are given by 


dp 

= 0 

dV z 

= 0 

dz 

dz 


dT 

= 0 

dB r 

= 0 

dz 

dz 

dV r 

= 0 

dB e 

= 0 

dz 

dz 


dVg 

dz 

= 0 

dB z 

dz 

= 0 


For the insulated boundary segment between ED we assign the boundary conditions 


Vr = V0—V z — O no slip boundary condition 


T = Tq = constant 



B r = B e = B z = 0 


For the uninsulated boundary segments FE and DC we assign the conditions 

V r = Vff = V z = 0 no slip boundary condition 
T = Tq = constant 




dBr 


dB z 


and simultaneously 


Jq — 0 which implies ^ 

d ±M = Q and ^ = 0 
dr dz 
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These later boundary conditions upon B insures that the electric field E satisfies the condition 
E ■ t = 0 everywhere on the nozzle boundary, where t represents a unit tangent vector to an 
arbitrary point on the nozzle boundary. 

Initial conditions are assigned in order to avoid large initial transients in the numerical solution 
because large changes can lead to numerical instabilities of the system of partial differential 
equations. We therefore assign the following initial conditions at all interior grid points. 

T = To = constant 

V r = V e =0 

V z = Vq — constant 

B r , Bq and B z are assigned values such that V • B = 0 
everywhere in the solution domain. 


TRANSFORMATION OF COORDINATES 


We make the change of variables 

z r 

x — y r j 7 

b f( z ) 

so that the domain of the nozzle 0 < r < /(z), 0 < z < b transforms to the computation domain 
0 < x < 1 and 0 < y < 1. The scalar form of the governing equations can then be written in the 

following forms. 

Continuity 

dp 1 djpVr) pVr 1 djpVz) vf'(z) djpVz) _ 0 

dt + /(z) dy + yf(z ) + b dx f(z) dy 

Momentum Equations 


dVr V^dVr_ fldVr _ vf'jz) 9Vy \ _ Yl 
dt + f(z) dy +Vz \bdx f(z ) dy J yf(z) 

dV 6 Vr dVp v (ldV 9 yf'(z)dV 9 \ VyVe 
dt + f(z) dy + z \b dx f(z) dy J yf(z) 


1 

P 




i=\ 


1 

P 


i=i 


dV z Vr dV z fl dVz 
dt + f(z) dy + 2 \b dx 


yf‘(z) dVz 

f(z) dy 
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Energy Equation 

dC p ( dT Vr dT fldT yf'(z) dT\\ 
p( c p+ t -qjT ) + /( z ) + z \bdx f(z) dy)) 


K t 


1 d 2 T 2yf '(z) d 2 T _ 

ft 2 dx 2 6/(z) ctedy 


, ( f'W 

K/( 


(f l_o(^ 

X) Z \f(x) 

2 (nz)\ 2 &T 1 _$T 1 &A+Y 4 . 

+y V M ) dy 2 + yf 2 ( x ) Q y / 2 ( x ) d y 2 ) h\ 

Electromagnetic Field Equations 


e r -component 


dT 

dy 


dBr f 1 d 2 b r 2 y£{z)^Br (yf^W 

~~df~ °\6 2 dx 2 bf(z) dxdy V f(z) J Q y 2 


+ 


f'(z) dB z 1 d 2 B z 


l_A 
01 { A 


f 2 (z) dy + bf{z) dxdy f 2 (z) dy 2 
T7 fldBr yf'(z)dB r \ , D (IdVz yJ\z)<W£ \ 

+ /(*) 9y) \l>9x f(z) dy ) 

/13B* _ (1 dVr yfXzj dV A 

Vr \bdx f(z) dy ) ‘\bdx f(z) dy J 

1 d 2 Bg 2yf'(z)d 2 Bg i ,yf'(z)^d 2 B l 

r ( y 


yf'(z) d 2 B ■, 


} 


-P 



b 2 dx 2 bf{z) dxdy V f(z) dy 2 


+ 


fldB z yf'(z) dlh \ / 1 dBg yf\z ) dbg \ 


\b dx f{z) dy J \b dx f{z) dy J 


f'{*) 

/(*) 


-2 


n*r 

/(*). 


dBe 

dy 


+ B, 


( f'(z)dB e , 1 d 2 B e y£{z)^B 0 1 (ldB, y/^)ggg\\ 

V f 2 (z) dy + bf(z) dxdy f(z) dy 2 yf{z)\bdx f(z) dy J J 


fldBr yf'(z) dB 
^ 1 b dx 


1 dB & B e 


-50 \ 

yf( z )J I 


f(z) dy J \f(z) dy yf{z) 

/I dT yf'{z) dT\( /l^_y/W!VV B (A 

^^\bdx f(z) dy J V \b dx f(z) dy J r \f(z) dy yf{z) 


dBe + Bq 
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e# -component 


dB e 

dt 




1 a 2 B # 2 yf’(z)d 2 B e yj'(z) 

bf(z) dxdy ' f(z) 


b 2 dx 2 


-y 


roo 

-2( f ' (z) ) 2 } 

/w 

\f(z)J J 


dB e 1 d 2 B 9 


dy 2 

1 dBe 


dy 


f 2 (z) dy 2 yf 2 {z) dy y 2 f 2 ( z ) ) 


Be \ 


yfjz) dBg\ 
f(z) dy ) 


(\ dV z yf'(z)dV z \ , V^dB^ ^ B^dV^ _ Ve_dBr_ _ B^^ 

+ B 0 T“5T"" I + f(r\ /)*/ + f(z\ dv f(z ) dv f(z) dy 


/»{ B 


1 d 2 Br 2 yf'(z)d 2 Br ( yf\z) ) 2 d^Br 

bf(z) dxdy K f(z) ' dy 2 

1 d 2 B z yf'(z) d 2 B z 


b 2 dx 2 

,f_{z)_dBz^ 

f 2 (z) dy bf(z) dxdy f 2 (z) dy 2 
( 1 dB z yf'(z) dB z \ ( 1 dB r yf'( z ) dB r 


- y 


f/"W- 2 | 



KK*)J J 


dBr 

dy 


\b dx f(z) dy 

2Be f l dBe_ _ yf'(z) dBe \ 

~ yf( z ) \ b dx f(z) dy ) 

n . - d 2 B r 

+ B t 


1 dB 


b dx f{z) dy f(z) dy 




(- 


yf'(z) d 2 B T 


1 d 2 B; 


+ 


(r. 


f ( z ) dy 
f3'(T ) dT 


/'(*) dBr T 

f 2 {z) dy bf(z) dxdy f 2 {z) dy 2 f 2 (z) dy 2 

dB r \ (IdBr y f'(z)dB r 1 dB z 
b dx 


B, 


/M dy ( 


l dBe 
b dx 


f(z ) dy f(z) dy 

yf'( z ) dBe 


f( z ) dy J 


— B r 


■)} 

1 dBr 
b dx 


) 




dBe 


+ 


yfP)dT\ 

/O) dy J 

Be ' 


B : 


1 dBr 


[b dx 


y_m dBr 
f( z ) dy 


Vf’(z) 

f( z ) dy 

1 dB 


1 dB , 


f( z ) dy 


f( z ) dy 


r) 


f{z) dy yf{z) 
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e z -component 


dB z _ ( f'(z) dBr 1 (PBr yfjz) &Br 

dT ~ a V PP) dy b f ( z ) dxd v / 2 ( z ) d y 2 

1 d 2 B z 1 p dB r yfj^dBA _ _Ji_dB z \ 
_ / 2 (z) dy 2 + yf(z) [b dx f(z ) dy J yf 2 (z) dy J 

V z dBr Br dv z V r dB z _ B z dVr VzBr | VyB z 

~J{P)~di~ f(z ) dy + f(z) dy + f{z) dy yf(z) yf(z) 

( f'(z) dB e , 1 d 2 B e yf'{z)d 2 B e \ 

+ / *VH PJ dy + bf(z)dxdy f(z) dy* ) 

J_dB 2 _ (1 dBg_ _ yJjpdBj 
+ f(z) dy \b dx f(z) dy ) 

B r d 2 Be B r dB e 1 dB r ( 1 dB e Bp \ 

+ / 2 (z) dy 2 + yp(z) dy + f(z) dy \f(z) dy yf(z)J 

B_z ( 1 dBp _ yf'(z ) dBp \ Br dBp \ 

+ yf(z ) \b dx f(z) dy J yP{z) dy J 


jT)dT r /lgfr 

/(«) i z \b dx 


yf'{z)dB e \ 
f(z) dy J 


+ B r 


( \ ^Be , Be A 

\f(z) dy yf{z)J 


NUMERICAL SOLUTION 

We axe primarily interested in the steady state solutions and the time necessary to achieve 
steady state. The above system of coupled nonlinear partial differential equations are simplified 
by assuming symmetry with respect to the 9 variable. This enables us to set all derivatives 
with respect to 8 equal to zero. Additional assumptions regarding magnitudes of force terms and 
energy terms can be made by doing a dimensional analysis of the resulting system of equations. A 
grid generation technique is used to alter the solution domain to a rectangle. Then the equations 
and rectangular boundary can be scaled before any numerical solution techniques are applied. 
The system of equations are then solved numerically using ADI (Alternating Direct Implicit) 
techniques patterned after the Lax modification. 

GRID GENERATION 

Let r = f(z) denote the nozzle boundary for 0 < 2 < b and consider the mapping from 
the (r, z) real coordinates to the (x,y) computational coordinates given by the transformation 
equations 

1 ” 6 r _ r < 23 > 
r max f( z ) 
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where r m ax = f( z ) denotes the nozzle boundary which changes with position. This mapping 
converts the region 

D = {r, z | 0 < z < b, 0 < r < r max } 
to the region D' of computational coordinates given by 

D' = {x,y|0 < x < 1, 0 < y < 1 }• 

To handle large gradients in any of the independent variables, the computational x,y domain is 
partitioned into 6 regions as illustrated in the figure 2. 
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Figure 2. Computational coordinates 


The 6 regions are characterized by the selection of the Ar and Ay step sizes. In this way finer 
grids can be specified near the boundaries and nozzle regions where large gradients can occur. 
Observe that any partial differential equation of the form 


du 

dt 


+ ^(r + D,{r,z)^ + D 4 (r,z + D b (r,z + D 6 


d z u 


du 


du 


(24) 
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where u = u(r,z,t ) and D 6 = D 6 (r,z,t,u , . . .) can be converted to computational coordinates 
x, y by using the chain rule for partial derivatives. These changes can be represented 

du du dx du dy 

dr dx dr ' dy dr 

du du dx du dy 

dz dx dz dy dz 


and 


d 2 u du d 2 x dx 

dr 2, dx dr 2 dr 

du d 2 y dy 

+ ~ O l” ^ 


d 2 udx d 2 u dy 


dx 2 dr dxdy dr J 
d 2 u dx d 2 u dy r] 


dy dr 2 dr [dydx dr dy 2 dr 

d 2 udx d 2 u dy 

dx 2 dz + dxdy dz J 
d 2 u dx d 2 u dy) 


d 2 u du d 2 x dx 

dz 2, dx dz 2 dz 

du d 2 y dy 
+ d^d^ + d^ 
d 2 u du d 2 x dx 

drdz dx drdz dr 


+ 


dx drdz 
du d 2 y dy_ 

dy drdz dr 


dydx dz dy 2 dz J 
d 2 udx d 2 u dy 


dx 2 dz dxdy dz] 
d 2 u dx d 2 u dy 


[ dydx dz ' dy 2 dz 

Then the partial differential equation (24) can then be written in the x,y computational 
coordinates 


du 


D i 

du 

d 2 u 

~dt 

— 

~7T~ x rr 

ox 

dx 2 



D 2 

' du 

d 2 u 


+ 

~^~ x rz 

dx 

dx 2 



D 3 

du 

d 2 u 


+ 

~77~ x zz 

dx 

dx 2 



Da 

r du 

du 


+ 

Tx Xr + di Vr 


. .o du 

{Xr) + faFy Zryr + V + dxdy 


d 2 u d 2 u 2 

yrX r + Q oCl/r) 


dx r 


d 2 u du d 2 u 

XrXz + fady XrVz + d~y Vrz + fody 


du 

y r Xz + ’fajVry 2 


, ,2 ^ 2 U , 

<**> + d^y X ’ V ‘ + 9^“ + dxdy 


+ Dq 


d 2 u d 2 u \2 

y z x z + -^-oiVz) 


dy 2 


- 

du du 

+ Ds 

~*~ Xz + o-y* 
dx dy 


or 


dt 1 


d 2 u 

dx 2 


+ Dl 


d 2 u 

dxdy 


+ D 


* d 2 u 
dy 2 


du 


du 


+ d:^ + d;^ + d 


dx 


dy 


(25) 
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where 


D\ = D*(x,y ) = £>i(x r ) 2 + D 2 x r x z + D 3 {x z ) 2 
D 2 = D%(x, y ) = 2D\x r y r + D 2 (x r y z + Vrx z ) + 2 D 3 x z y z 
D$ = Di(y r ) 2 + D 2 y r y z + D 3 y 2 

D\ = D^(x, y) — D\x r r + D 2 x rz + D 3 x zz + D±x r + D$x z 
D% = Dl(x,y) = D\y rr + D 2 y rz + £>31/ zz + I>4j/r + D$y z 
D§ — D$(x,y) = DG{yf{x^,x,t,u,...) 
and Di = Di(r, z ) = £> t (y/(x), x) for i = 1, . . • , 5 and 


x r =0, y r = 1 //(*), = £, Vz = - yf'{z)/f(z ) 

X rr = x rz = x zz = y rr = Q, yrz = -f'( z )/(f( z )) 2 i yzz = ~V 


Af) 2 /Af) 
/(*) v /(*) 


2 


Then the above coefficients reduce to 


„ n f(z) 

D 2 (x,y) = 7T7-T -2yD 3 - 


D 3 (x,y) = 
D\(x,y) = 


b f( z ) 
D i 


bf(z ) 

/'(*) 


(/(*)) 

£>5 




2 n /7'MY 

v Itw j 


DU x > y ) = - g 2 (^( 1))2 - s '® 3 

£> 6 (*, y) = £^6 {yf( z )i x,t,u, . . .) 


f_ 

/( 


!M _,/A£>V'\ 
fM l/wl / 


+ ^W' W 


/(*) 


/M 


where 


Dj = D{(r,z) = Di(yf(z),x), for i — 1,2, 3, 4, 5. 


ADI NUMERICAL METHOD 

The following description of the ADI numerical method is for uniform Ax and Ay grids and 
of course has to be modified for unequal x and y spacing. In step 1 of the ADI numerical method 
the interior points to the region D 1 of the computational domain are labeled from left to right as 
illustrated in the figure 3. Assume that the system of partial differential equations to be solved 
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have all been normalized. Partition the segment from x = 0 to x = 1 into segments with spacing 
Ax = l/m 2 so that the ith node is iAx and the right boundary is m 2 Ax, with 0 < * < rn 2 . 
Similarly, partition the segment from y = 0 to y = 1 into segments with spacing Ay = 1/mi so 
that the jth node is jAy and the top boundary is mi Ay with 0 < j < mi, where i,j,m\ and 
m 2 are integers. The interior points to this grid are then labeled as illustrated in the figure 3. 
Let u n be associated with the ( i,j)th node point, where 

n = (m 2 - l)(mi - 1 - j) + i mi and m 2 are fixed. 

Conversely, given the u n point, we can determine its position i,j from the relations 

j = mi — 1 - Int[(n - l)/(m 2 - 1)] 
i = n - (m 2 - l)(mi - 1 - j) 

where Int[x] is the integer part of x . 


u 1 

u ni2 


U2 

u rri2 + 1 


U 3 U 4 

^ m 2~}“2 ^ m 2+3 


U YTl2 — 1 

u 2rri2 —2 


U„ 

u (m2-l)(mi — 1) 



Figure 3. Step 1 labeling of interior points to domain D 
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Letting 


u(iAx,jAy,nAt ) = u"j 

and dropping the star notation, we then replace all partial differential equations of the form of 
equation ( 31 ) by difference equations having the form 


t^f 1 _ 

«.J Id 

At 


= D i 


+ D 2 


+ Dj, 


+ D§ 


u^-2u?r+uv 


"+ 1 + u n+ .' . 


(Ax)- - 


“"■M.i+l -“"+U-1 


4 AxAy 




(Ay)" 


) 


+ Da 


# n+l _,, n + 1 .' 

h+l,i t-l,j 

2Ax 


U i,j + 1 


i -"k-A 

2Aj/ / 


+ Dq 


which can be rearranged to the form 




2£>iAA ^ n+ i 
(Ax) 2 / 


u * ' — 


\(Ax) 
/£> 3 A< 
' \(Ay) 2 


-D4AA n+ i 


2 + 2Ax 


A 


t+i,i 


: + 


L>i At \ n+ i 


£> 5 a*\ 

2Ay ; 


/ £>4 At 
V 2 Ax 
(£3 At _ £ 5 AA u „ 


(Ax) 2 ; 


ij 


( 2D*At\ „ 

= (' - j “<■>• + l(A^ + Wj Uij+1 V(A!/) 2 2Ay ) 

+ (<n,,- + i -< + w-. -“?- 1J+ i + “?-u-i) +^ A ' 




Evaluating this equation at each of the interior node points gives rise to a system of 
(m 2 - l)(mi - 1) implicit linear equations which are then solved by row reduction methods. 

The second step of the ADI method relabels the interior points of the computational domain 
from top to bottom as illustrated in the figure 4 . 

For this labeling we can let u n denote the point associated with the ( i,j)th node where 


n = (mj — l)(i — 1) + mi — j. 


Conversely, given u n we can solve for i and j from the relations 

i = 1 + Int[(n - l)/(m\ — 1)] 
j = (mi - l)(i - 1) + mi - n 
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Ul 

u m\ 

u 2 

w mj +1 

U3 

U mi +2 

U4 

u mi-j-3 

mi —1 

u 2mi —2 



Figure 4. Step 2 labeling of interior points to domain D’ 


For step 2, all the partial differential equations of the form of equation (31) arc 
difference equations having the form 


At = D ' l (SF 


r ^+u+i ~ u h-i,j-i ~ u ^-i j+i + a 
+ L>2 V 4AxAy ) 

/«?., • - u? 

D, { 2A^ 


u »,j+i _iA _iizl 

(Ay) 5 

n+1 _ n+1 

+ AH lJ+ L ^ l+At 


+ l ' 


2Ay 


replaced by the 
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which can be rearranged to the form 


(■* 

0 


(Ay) 2 / ' 3 


, 2DiAt\ 
- 1 1 " (Ax ) 2 ) 


3 ~ 


\(A y ) 2 

/DiAt 

V(Ax) 2 


2A y ) 

d 4 aa 

+ 1atJ 


+ 


P 2 At 

4AxAy 


(<- 


+i,j+i “ u «+i j-i u *-i,i+i 


n+1 , /'^At _ AjAA 

•J+ 1 + \ 2Ay (Ay) 2 ) 
n /Pi At D 4 At 

u i+i,j + Y( Ax) 2 2 Ax 

+ 


u 


n+l 


u 


n 


Applying this difference equation to each interior node point results in an implicit system of 
( mi - i)( m2 - 1) simultaneous linear equations which must be solved for the values of u at the 

node points. 

One can see that the finer the interior grid there results a larger system of linear equations 
to be solved. Also the results from the ADI numerical method are more accurate on the even 
numbered time steps. Additional complications results when employing the unequal step size 
approximations illustrated in the figure 5. The unequal step sizes are necessary to handle large 
gradients occurring in any of the dependent variables. The computational region is therefore 
divided into 6 regions as illustrated in the figure 5. The density of node points can be changed 
in each region by selecting different step sizes in the computational coordinates. 



Figure 5. Unequal x and y spacing. 
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In the case of unequal grid spacing we employ the difference approximations. 


du 1 



dx Aa:i + Ax 2 

du 1 



dy Ayi + At/2 

d 2 u 2 


dx 2 Axi + Ax 2 

d 2 u 2 


L |^(u(x + Axi , y) - u(x, y)) - ^(u(x “ Ax 2, y) - “(*> y)) 

^-(u(x,y + Ayi) - u(x,y)) - ^-(u(x,y - Ay 2 ) - «(x,y)) 
Ayi ^J/2 

~ u(x + Axi,y) - u(x, y) u(x - A x 2 ,y) - ^( x >y)) 

Axi Ax 2 

u(x, y + Ayi) - u(x, y) u(x,y - Ay2) -u(x,y)) 

Ayi A y2 


dy 2 Ayi + Ay2 
a 2 tx _ ujx + Axi , y + Ayi ) - u(x + Ax t , y - Ay 2 ) 
dxdy ~ (Axi + Ax 2 )(Ayi + Ay2) 

u( x - A x 2 ,y - Ay2) - u(a - Ax 2 ,y + Ayi) 
(Axi + Ax2)(Ayi + Ay2) 


+ 


SPECIAL CASE-ELECTRIC FIELD IN VACUUM 


In a vacuum we solve 


« d 2 <j> 1 d<f> d 2 <f> _ 

v ^ = ^ + ;¥ + a ? =0 


over 


the domain 0 < r < f(z) and 0 < 2 < 1. Using the transformation equations 


(26) 


Z = X 


r = yf(z) 


where 


is 


f(z) = .2 + x tan(87r/180) 

used to describe a straight line nozzle boundary. The equation (26) then transforms to 

d 2 <j> . , sd<f> 


d 2 * ^ , d 2 <f> 

+ a ( x ) y)‘ 


0 X 2 1 v ’ dxdy 
over the domain 0 <x<l, 0 <y<l where 

n*) 


+ 6(x,y)^2 +c(x,y)^ = ° 


(27) 


a(x,y) = -2y 


/(*) 


. 1 . ( /'(*)' 2 
K x ,y) = 72 + » y ~ 


P(z) 


/(*) 


, , 1 
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VOLTAGE 


F/4t/?£(,, 


POTENTIAL FIELD IN VACUUM 







flmef. 

MAGNITUDE OF ELECTRIC FIELD IN VACUUM 



LTAGE 







Using the notation U{ j = t/(zAx,j Ay) and the difference approximations 


Xlxx 


U X y 

u yy 

Uy 


uj+rj ~ 2 *+,j + M «-i J 

h? 

u+t-lj+l ~ ««•+!, j-1 ~ + «»-l j-1 

4/i 2 

u »,j+ 1 ~ 

h 2 


. u «+i,i ^izlA 
2 h 


and defining 

d( x ,y) = + K x iJ/)) 

the equation (26) reduces to the difference equation 

«t,i = ^77 | t+1 ’ J '^2 * ~ + 4^ ( u i+l,i+ 1 “ u *+l.i-l ~ u *-lj'+l + u i-l,j-l) 


U IJ 

h 2 


^ (“»,j+l + + 2h ( U,+1 ’ J U *- 1 J‘)} 


This difference equation is subject to the boundary conditions 

u ijmax = assigned potential value 
U 0J == U 2J 

^i,0 = ^i,2 

which represent zero derivative boundary conditions along the other three sides. The figures 6,7,8 
and 9 illustrate the potential function for two different nozzle configurations where the cathode 
is assigned a value of -500 volts and the anode(s) is assigned a value of +500 volts. 


FLOW AND HEAT TRANSFER THROUGH A POROUS MEDIA 


In the figure 10 a porous material is heated with 40kw of power from a solar simulator. We 
assume that the solid porous material is heated to a uniform temperature T s and that a gas flows 
through the porous material and is heated. 
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Q rad cooling 



Figure 10. Heat transfer through porous material. 



In the following discussions we use the notations: 

<f> = porosity 0 < <f> < 1 
T g ,T s = Temperature of gas and solid (K) 
p g , p s = Density of gas and solid ( gm/cm ) 

K g , K s = Thermal conductivity of gas and solid ( cal/s cm 2 K/ cm) 

Ke = (1 — (j>)K s + <f>K g = Effective thermal conductivity 
u = Velocity of gas (cm/ sec) 
a g = KgjC vg pg = Thermal diffusivity of gas (cm 2 /sec) 

L = Thickness of porous material (cm) 
uL 

Pe -- — = Peclet number ( dimensionless ) 

a 9 

r = Radius of disk (cm) 
h = Heat transfer coefficient (cal/s cm 2 K ) 

Qo = Input power (Kw) 

A = Surface area of disk (cm 2 ) 

Cpg,C ps = Specific heat of gas and solid (cal/gmK) 

R s = Ratio of surface area to volume of porous media (cm 2 /cm 3 ) 

U = T s /T g o = Dimensionless temperature ratio 
V = T g /T g o = Dimensionless temperature ratio 
r = — = Dimensionless time 

L 

x 

X = — = Dimensionless distance 

L 

Following reference 1, the basic equations describing the heat transfer to a gas moving through 
a porous media are given by 

^{w + ^) =K °i£- hR{T °- T ’ ) (28) 

Pa C p , ^ = A-,0- hR(T, - T,) (29) 


for 0 < x < L and t > 0. 
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The equations (1) and (2) are subject to the boundary conditions 


dT s | = 239 Qq 

: dz ' 2=0 ” A 


— ecr(Tg — Tgo) — pgC pg u(T s T g o) 


-K s ^ = h(T s -T w ) (31) 

or 

where all terms have been scaled to the units of cal j s cm and 

a = 5.67(10" 12 )(.239) cal/s cm 2 K 4 

is the Stephan-Boltzman constant, 6 is the emissivity. In addition we have the boundary 


conditions 


PgCpg U Qz L=0 ~~ hR(T S Tg) 


“•I £ =0 and flu =0 

dz 1 Z=L dz ' Z ~ L 


We further assume that the initial conditions are 


T s — T g — TgQ. 


Introducing the dimensionless variables U,V,X,t, the above equations can be written 


du _ . (\ d 2 u i du j_3Hr\ _ ... _ 

dr ~ A ° \r 2 dR 2 r 2 R OR + L 2 HZ 2 ) 


dV dV f 1 d 2 V 1 dV 1 d 2 V 
1fr + -dZ~ M + r 2 RdR + L 2 dz 2 


^ -B l( y-C) 


k. 


<v- 1.) 


where 


4 x 2 D h Rd? 

Al ~ Pe Bl Pe K a 


A 0 = 


A',X 2 


£o = 


ZiiLX 2 


= K DX -TTKg ~°-a gPs C ps Pe a gPs C ps Pe 

The equations (35) and (36) are subject to the boundary conditions 


= ^3 - B 3 (U 4 - I/ 4 ) - C 3 (l/ - V) 

%\ Z= o = B t{ U-V) 

_ o — o 

^7IZ=1“ U ’ 1^ = 1 U 
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where 


^3 = 


239 Q(r 0 R)L 


AK e T g0 ’ 

The initial conditions are 


#3 = 


™t* 0 l 

Ke ' 


KgPe 

c * = n<7 


B 4 = 


hRJ? 

a g Pe PgCpg 


U(0,t) = l and V(0,t) = l 


(41) 


(42) 


We divide the intervals 0<Z<land0<J?<l into N parts with step sizes 
AZ = A R = 1 /N. We desire to represent the temperature of the gas and solid at the positions 
R = iAR and Z = jAZ for the time r = nAr which is based upon the given temperatures 
at time r. Let U (iAR JAZ, nAr) = and V{iAR,jAZ, nAr) = Vg, and use difference 
approximations to write the above equations as difference equations. We use the ADI (Alternating 
Direction Implicit) method to solve the above system of coupled partial differential equations. 
Material Properties 


From reference 15, we obtained the following empirical data for Hafnium carbide. 


Temperature 
deg K 

Thermal Conductivity of Hfc 
W/cm K 

560 

0.09 

800 

0.12 

1100 

0.13 

2000 

0.15 

2500 

0.25 

3000 

0.29 


The best fit second degree polynomial to the above data is given by 

K,(T) = .0361887 + 1.03093(10 -4 ) T - 1.2077(10 -8 ) T 2 
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Temperature 
deg K 

Specific Heat of Hfc 
cal/gK 

60 

0.011 

95 

0.020 

300 

0.045 

600 

0.065 

1000 

0.065 

3000 

0.065 


Table lookup will be used to fit the Specific Heat data. 


Temperature 
deg K 

Emittance 

300 

1.00 

1100 

0.98 

1900 

0.90 

2500 

0.70 

2900 

0.62 


The above data is represented by the approximating function 

e(T) = .8 - .2 tanh((T - 2100)/1000). 

We use the above data to construct empirical relations to represent the constants in the above 
system of coupled partial differential equations.The Appendix A contains graphical output from 
the computer analysis of the heat transfer in a porous media. 
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APPENDIX A 

GRAPHICAL REPRESENTATION OF RESULTS 


Solution of the coupled equations describing flow of gas through a porous media. The 
input power in Kw is assumed to be a cosine curve of the form 

ttR 

Q(R) = 40. cos(— ). 

Solution of coupled equations is by ADI (Alternating Direction Implicit) technique. 
All dimensions have been normalized. The radial and axial directions range from 0 to 1 . 
The temperatures of the gas (V) and solid (17) have been normalized by the equations 

U = TS/TGO V = TG/TGO 


where TGO = 300 K. 

The figures Al throug A16 illustrate the temperature change for the gas and solid as 
a function of the normalized time 

t u 

t = t 

where t is real time, u is velocity, and L is length in the axial direction. 
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Figure Al. Normalized Gas Temperature at r — 1.0 
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Figure A2. Normalized Gas Temperature at r = 
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Figure A3. Normalized Gas Temperature at t — 3.0 
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Figure A4. Normalized Gas Temperature at r = 4.0 
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Figure A5. Normalized Gas Temperature at r = 5.0 
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Figure A6. Normalized Gas Temperature at r = 10.0 
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Figure A7. Normalized Gas Temperature at r — 20.0 
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Figure A8. Normalized Gas Temperature at r = 30.0 
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Figure A9. Normalized Solid Temperature at r 1.0 
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Figure A10. Normalized Solid Temperature at r = 2.0 
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Figure All. Normalized Solid Temperature at r = 3.0 
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Figure A12. Normalized Solid Temperature at r — 4.0 
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Figure A13. Normalized Solid Temperature at r = 5.0 
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Figure A14. Normalized Solid Temperature at r — 10.0 
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Figure A15. Normalized Solid Temperature at r = 20.0 
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Figure A16. Normalized Solid Temperature at r 
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